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O ■ Abstract 

■ We have studied the statics and dynamics of flux hnes in a model for 

a 

I ' YBa2Cu307_5 , using both Monte Carlo simulations and Langevin dynam- 

^ ' ics. The lines are assumed to be flexible but unbroken in both the solid and 

O 

liquid states. For a clean system, both approaches yield the same melting 
^ ■ curve, which is found to be weakly first order with a heat of fusion of about 

0.02kBTm per vortex pancake at a field of 50kG. The time averaged magnetic 
field distribution experienced by a fixed spin is found to undergo a qualita- 
tive change at freezing, in agreement with NMR and /xSR experiments. The 
calculations yields, not only the field distribution in both phases, but also an 
estimate of the measurement time needed to distinguish these distributions: 
we estimate this time as > 0.5/xsec. The magnetization relaxation time in a 
clean sample slows dramatically as the temperature approaches the mean-field 
upper critical field line Hc2{T) from below. Melting in the clean system is 
accompanied by a proliferation of free disclinations and a simultaneous dis- 



appearance of hexatic order. Just below melting, the defects show a clear 
magnetic-field-dependent 3D-2D crossover from long disclination lines paral- 
lel to the c-axis at low fields, to 2D "pancake" disclinations at higher fields. 
Strong point pins produce an energy varying logarithmically with time. This 
In t dependence results from slow annealing out of disclinations in disordered 
samples. Even without pins, the model gives subdiffusive motion of individual 
pancakes in the dense liquid phase, with mean-square displacement propor- 
tional to t^/^ rather than to t as in ordinary diffusion. The calculated melting 
curve, and many dynamical features, agree well with experiment. 
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I. INTRODUCTION 



There have been numerous debates about the superconducting-to-normal transition in a 
magnetic field. Abrikosov's classic mean field theory predicts a rigid vortex lattice persisting 
to the mean- field transition line, Hc2{T). More recent work suggests that the mean- field 
Abrikosov transition is altered by fluctuations in both two and three dimensions (2D and 
3D). In the 2D case, numerical evidence suggests a first order melting transition from a 
vortex solid to a vortex liquid , followed at higher temperatures by a gradual crossover 
to a normal state. In 3D, the first suggestion of a first-order fiuctuation-induced transition is 
due to Brezin et a/ @] . Their argument is based on Landau- level expansion of the Ginzburg- 
Landau functional for a conventional superconductor. 

Much recent work on 3D fiux-lattice melting transition has been stimulated by the be- 
havior of high-Tc materials . Several experiments have convincingly demonstrated that 
the transition there is first-order in a sufficiently clean system at sufficiently low fields . 
On the theoretical side, numerical studies of melting have been based on model pairwise in- 
teractions ||T0Hl5[, the frustrated XY model [p!^-[T8[|, and an expansion of the free energy in 
lowest Landau levels (LLL) 0J^, among other approaches. Some of these calculations PJT7[] 
also indicate that the melting transition in the pure 3D system is first order, as suggested 
by experiment. 

However, numerous issues remain unresolved. One such issue is the effects of disorder 
on the first-order transition. Disorder is widely expected to modify the first order transition 
to either a vortex-glass [|l9l or Bose-glass transition [^, depending on the type of disorder. 
Modeling such disorder is difficult because of problems associated with slow relaxation. 



Much of the modeling has therefore been carried out within the frustrated XY model [21 



a model which introduces an artificial pinning by a fictitious lattice. Other unresolved 
issues center on the dynamics of the solid and liquid vortex system, for which only a very 
limited number of calculations have been carried out. Dynamical calculations are obviously 
necessary to understand many measurable properties of the high-Tc materials, such as the IV 
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characteristics, voltage noise spectra, NMR, and ^SR. A common approach is to model the 
high-Tc material as a network of Josephson junctions (a natural dynamical generalization 
of the frustrated XY model) |^2|]. But as in the static case, this model also suffers from the 
problem of fictitious pinning, though it is reasonably tractable numerically. 

A few dynamical calculations have been carried out using a time-dependent Ginzburg- 
Landau (TDGL) model within a vortex representation. For example, Enomoto and cowork- 



ers have studied various aspects of flux lattice melting using this approach. They focused 
on the effects of pin density on the irreversibility line, using as a melting criterion the on- 
set of flux line diffusion. They also considered the transport properties of pinned systems, 
but studied only configurations consisting of either a single flux line or a two dimensional 
lattice with point disorder. Reefman and Brom applied a similar technique to a model for 
a single layer of BiSr2Ca2Cu208 , neglecting the Josephson coupling between layers and fo- 



cusing on the NMR properties Most recently, Probert and Rae |]T5| have performed a 
Langevin simulation for a model YBa2Cu307_5 system, both with and without pins. Their 
results shed some light on the dynamical distinctions between the irreversibility line (found 
in disordered systems) and the thermodynamic melting line (characteristic of a clean sys- 
tem). Their calculations, however, assume rigid vortex lines, thereby leaving out some of 
the most characteristic three-dimensional (3D) behavior associated with flexible lines. All 
these results suggest that this approach may be a useful and quite realistic way to treat the 
dynamics of flux lines in high-T^ materials.. 

In this paper, we present a numerical study of both the statics and dynamics of a "layered 
London model" for a three dimensional flux lattice in YBa2Cu307_5 , using a combination of 
Monte Carlo (MC) simulation and Langevin dynamics (LD) within a vortex representation. 
Our work extends earlier studies in a number of ways. For example, we determine the melting 
line not only by the motion of individual vortex lines, but also by changes in equilibrium 
quantities such as the vortex structure factor and a hexatic order parameter. We also 
determine the conditions under which the vortex lines maintain their integrity, even in the 
liquid state. Perhaps of greatest interest, we find a novel 3D-2D crossover in the structure 



of the vortex solid just below melting, at which the characteristic topological defects change 
from long disclinations parallel to the c-axis to short disclination "pancakes." This crossover 
may possibly be connected with some recent experimental work, as discussed further below. 

The remainder of this paper is organized as follows. The next section describes the model 
and discusses our choice of parameters suitable for YBa2Cu307_5 . The following section 
presents our numerical results as obtained by both Monte Carlo and Langevin simulations. 
A brief discussion follows in the final section. 



II. MODEL 



A. Model Classical Action 



In our model, interest is confined to fluctuations in the phase degrees of freedom of the 
superconducting order parameter ip. The amplitude {ipl is assumed not to fluctuate, but 
instead takes the value dictated by minimizing the Ginzburg-Landau free energy at the 
given temperature T and magnetic induction B. This resulting is related to an effective 
in-plane penetration depth p3| by 



m c 



A^.(0) 



where T'c(O) is the mean- field transition temperature at zero magnetic field, Aa;,(0) is the 
in-plane penetration depth at zero temperature, and Bc2(T) is the mean-field upper critical 
field line. is normalized so that m* is twice the electron mass and e* = 2e. 

The model consists of Nz parallel superconducting layers a distance d apart. Each layer 
contains N^^ two-dimensional "vortex pancakes" (i. e., 2D vortices) described by transverse 
position coordinates rj^^. The vortex density is assumed fixed at = l/a% = B/(f)o, 
where 0o = hc/2e is the flux quantum. Such pancakes in different layers interact via both 



magnetic and Josephson interactions [pl| -|27[]. Here, we simply assume that these interactions 
combine to produce flexible but unbreakable vortex lines - that is, each pancake is always 
associated with two specific pancakes in the adjacent layers. The justification and possible 



limitations of this assumption are discussed below. The interlayer coupling strength is 
characterized by a single variable 7 = ^a6(0)/^c(0), where ^ab{0) and ^c(O) are the zero- 
temperature superconducting coherence lengths in the ab-plane and c direction. 7 has 
associated with it a length scale Vg = '-yd. The layered structure becomes important for 
lengths shorter than [EH . 



Following |Tl|, we write down the Hamiltonian for the system as 



i k 



\^i,k '^j,k\ 
2^n 



Here the in-plane repulsive interaction takes the form 



U{x) 



while the interlayer interaction is taken as 



(2) 



(3) 



V{x) = cj{x - 1) {x > 1); 



Cj{X 



ix < 1] 



(4) 



with 



87r3A„,(T,i?)2 



1 + ln 



Aafe(O) 

d 



(5) 



In order to reduce finite size effects, we employ periodic boundary conditions in all directions. 
Because of these, the effective in-plane interaction becomes Kq{x), which represents the 



summation of the modified Bessel function Kq{x) over image vortices ||29[| 



B. Langevin Dynamics 



To probe real time dynamics, one can also run LD simulations on the same model, 
assuming that the vortices are subject to an overdamped dynamics. Then the equation of 
motion for a vortex pancake can be written 



i,k 



(6) 



The first term on the right-hand side is the Brownian force due to thermal noise. The noise 
is assumed to be Gaussian-distributed white noise with correlation functions 

(/SW) = 0' (7) 

((fj(t) ■ h^) iq,,it') . hp)) = 2kBT'l5,A,p5,y5{t - t'), (8) 

where ha is a unit vector in the a direction, a = x,y. The second term on the right-hand 
side of eq. (|) is the force due to the other pancakes; it is obtained as a negative spatial 
gradient of the vortex- vortex interaction term, as written down in eqs. The third 

term (not studied numerically in the present paper) is the Lorentz force due to an apphed 
current. 

The last term in describes the force due to the random pinning potential. The pins 
are modeled as uniformly cylindrical regions of radius [taken to be 2^ab(0) throughout this 
work] . The pinning energy of the vortex pancake is expected to be determined by the fraction 
of the core area within the pinning well. We can achieve this dependence by assuming a 
pinning energy per pancake Up(T,B) = ap(i0g/[167r^A^^(T, _B)]. Thus, the strength of a 
single pin is controlled by the dimensionless parameter ap. The effectiveness of the pins is, 
of course also influenced by their areal density Up, or the "equivalent field" Bp = (poUp. 

For simplicity, the force due to the l^'^ point pin is assumed to be directed radially inward 
towards its center (at R;). For Vp > C,ab(T) it is given by 

fF{^i,k) = -^Wl^ if rp-a(T)<|r,fc-Rz|<r, + a(T), 

= otherwise. (9) 

For Vp < ^afc(T), it takes the form 

rp + Ub{T) Hafe(T)^ 
= otherwise. (10) 

This choice includes in the simplest manner the fact that the vortex core area grows with 
increasing T while the disorder is temperature-independent. 



C. Numerical Approach and Choice of Parameters 



To obtain the thermodynamics via MC, we use the standard Metropohs algorithm with 
variable step sizes, as discussed in |rT]]. Typically, we equilibrate over 2 x 10"^ MC steps and 
evaluate the thermodynamic averages over an additional 2 x 10^ — 10^ steps. 

For both MC and Langevin calculations, we use look-up tables for both the potential and 
the forces, as well as a scheme for interpolating between the points in the table. The time 
iteration is carried out using a second-order Runge-Kutta algorithm in time steps of A ■ to 
where to = ^j^^r^ and = g^2^'^°(o-)2 • The choice of A is dictated by the dominant forces 
in our model [eq.@], which in this paper are the vortex- vortex interactions; in general, 
A < (9(10). For optimum convergence, we have allowed A to depend somewhat on T,B, 
and J, since different components of the force may dominate at different values of these 
parameters. 

Finally, we briefly discuss our choice of parameters. For Aab(O), we use 1000 A for 
YBa2Cu307_5 single crystal. This choice is close to the experimentally determined value, 
and also places the simulated melting curve in close agreement with the experimental data 
of IP). For the remaining parameters, we use the values appropriate for YBa2Cu307_5 : 
K = 87.5, d = 11. lA, 7 = 5, Tc(0) = 93K, and dHc2{T)/dT = -1.8 x lO^Oe/K, where 
K, = Aab(0)/^afe(0) is the Ginzburg-Landau parameter. 



D. Calculated Quantities 

Before discussing our numerical results, we first define a few important physical quanti- 
ties. 6R{t) is the transverse root-mean-square (rms) displacement of a pancake vortex from 
its average position, i. e. 

mt) ^ - l^T.{i^^,>^- < >*)'>!^'| ■ (11) 

Here < ... >t denotes an average over time t, and Ntot = NzN^ is the total number of 
pancakes. Now in a finite system, as in our simulation, the collection of vortex lines tends to 
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drift as a whole even in the sohd phase. This behavior, seen in both sohd and hquid phases, 
is strictly a finite-size effect and has no relation to any measurable quantities. We therefore 
subtract out the drift by using r*^(t) = Ti^kit) — Hcmit) instead of rj^fc in Eq- (O), ^cmit) 
being the instantaneous center-of-mass coordinate of the entire lattice. To monitor lateral 
fluctuations, we also calculate the "wandering length" It deflned by 

^T = JJ^T.{\^i,k-ri,k+i\^) (12) 
In addition, we compute the density- density correlation function 

C(r, z) =< p^(r, z)p^{0, 0) >t^oo, (13) 
and its partial Fourier transform 

S(q..)=/dVC(r..)expe,q.r), (14) 
where Pt,(r, z) is the local vortex number density in each plane z and q = {q^, Qy). 



III. NUMERICAL RESULTS 
A. Location of Melting Curve; Heat of Fusion 

To locate the melting point for a given B by MC simulation, we flrst make a quick sweep 
over a wide range of temperatures (~ 20K), using 2 x 10^ MC steps for each T in steps of 
AT ~ 0.5 — IK. We interpret a discontinuous jump in 6R{t) as well as the vanishing of the 
intensity of the Bragg peak ^(q, 0) at q = G2 (where G2 is a 2D reciprocal lattice vector of 
the triangular lattice) as signatures of melting. We then repeat more careful sweeps over a 
narrower temperature region with up to 10^ MC steps. 

If the temperature is cycled through using an interval of 10^ — 3 x 10^ MC steps for 
each temperature in increments AT/T^ = 0.0024, we observe hysteresis in most monitored 
quantities. The width of the loop is typically ~ O.OIST^. Hysteresis is most pronounced for 
6R{t) and for the disclination density (deflned below), but is also quite conspicuous for the 
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hexatic order parameter(also defined below) in the same temperature range. From the size of 
the jump in the total internal energy seen at melting, we estimate the latent heat per vortex 
pancake to be about 0.034 ± O.Ol/^^Tm. (A possibly more precise estimate is given below.) 
The melting transition is calculated to occur at Tm{B)/Tc2{B) = 0.87 ± 0.02, 0.93 ± 0.006, 
and 0.93 ± 0.02, for B = 90, 50, and 10 kG respectively, and lattices with eight layers, 
in reasonable accord with experiment [H. For lattices with 32 layers and Xab = 1400 A, 
we obtain T^{B)/Tc2{B) = 0.86 ± 0.02,0.92 ± 0.02, and 0.96 ± 0.02 at the same fields. 
[Here Tc2{B) is the mean-field transition temperature at field B.] Since our assumed T- 
dependence of Xab{T,B) is likely to become less accurate as T — > Tc(0), we may expect 
increasing deviations from experiment at lower fields, as indeed we find in our calculations. 
Note that, for the higher fields, Tfn/Tc2{B) depends little on the either the lattice aspect 
ratio or the value of Xab- This behavior is consistent with the dimensional crossover discussed 
below. 

The transverse displacement 6R{t), as expected, shows characteristically different be- 
havior in the solid and liquid phases. In the liquid, SR{t) increases with increasing t (see 
below), while in the solid phase, it saturates after a short transient period. The inset to Fig. 
1 shows the behavior of 6R{t) across the melting transition for B = 50kG, as calculated by 
MC. Also shown are the Langevin results, which agree very well with MC predictions. This 
confirms that the two indeed give, as expected, very similar predictions for thermodynamic 
quantities. From the results in the solid phase, we can read off the Lindemann number 
cl = 6R{Tfn, B) = 0.18, at B= 50kG. We have not, however, confirmed that cl is constant 
along the melting curve, as would be expected if Lindemann's law is really valid. 

Another way to display the first-order melting process is also shown in Fig. 1. The MC 
energy at any given temperature fluctuates with MC time. The two curves in the main part 
of Fig. 1 show the distributions of total internal energy within two different time windows at 
the melting temperature [(a) and (b) in Figs. 1 and 2]. On either side of T^, this distribution 
tends to be independent of the initial time of the window. Precisely at T^, however, the 
system slowly oscillates between a liquid and a solid phase with a correlation time of about 
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5 X 10^ MC time steps (for this system size). The oscillation causes the two distributions to 
differ: the two distributions shown are the extreme cases that we have found in the length 
of time studied. As shown in the insets at the bottom of Fig. 2, the density correlation 
functions in the "low-energy" and "high-energy" windows do indeed show solid-like and 
liquid-like characteristics. We believe that the weak hexagonal symmetry seemingly present 
in the "liquid" phase actually results from the presence of a few solid configurations in the 
window, and possibly also from the finite size of the sample. 

/^From the average energy difference between the two phases as read from Fig. 2, and 
also from the distance between the peaks in Fig. 1, we estimate the heat of melting to 
be ~ 0.02kBTm per vortex pancake. This estimate for the heat of fusion is in agreement 
with LLL calculation of 0, as well as with experimental estimates at high fields. It is 
significantly smaller, however, than Hetzel et al U^'s value of O.SkBTm per vortex pancake, 
obtained using a frustrated stacked triangular XY model. The reasons for this difference are 
a matter of speculation. One possibility is that in both our simulations and the LLL model, 
which sets the overall energy scale, decreases with increasing temperature, whereas in 
the XY simulation, the coupling strength J is temperature- independent. More plausibly, 
our simulations are carried out at high fields while the XY simulations are more appropriate 
at low fields. Yet a third (unlikely) possibility is that our simulations allow for anisotropy 
in the superconducting properties, whereas those of Hetzel et al do not. 

We find that most of the energy discontinuity at melting comes from changes in the 
interaction energy between different pancakes in the same layer, which has a rather clear 
jump at melting. By contrast, the interlayer coupling energy has large fluctuations in both 
the solid and liquid phases but no clear jump. These fluctuations tend to mask the step-like 
change of the in-plane component. 
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B. Topological Defects in Clean YBa2Cu307_5 



We have carried out a search for topological defects above and below melting in our 
model for YBa2Cu307_5 . Basic building blocks for various kinds of topological defects 
are disclinations and dislocations. Both kinds of defects may be identified by carrying 
out a Delaunay triangulation |^| for individual layers of each sampled configuration. To 
understand this procedure, we show in Fig. |^ a snapshot of a typical melted configuration 
in a single layer which shows both bond lines, as identified by the Delaunay procedure, 
and topological defects. The vortex pancakes are located at the vertices of the triangles. 
They are marked by black and grey dots if their number n{i) of in-plane neighbors is seven 
(positive disclination) or five (negative disclination) rather than the six expected for a perfect 
triangular lattice. We see characteristic examples of an isolated dislocation [i.e. a pair of 
disclinations, marked by (a)], an isolated disclination (b), and a bound pair of dislocations 
(c). In the present work, we arbitrarily call a pair of disclinations "bound" if and only if 
they reside on neighboring pancakes as in (a) and (c). 

For a more quantitative analysis of the defect configurations, we assign "charges" q{i, z) = 
n{i, 2;) — 6 to the iz^^ pancake, where n{i, z) is the number of in-plane neighbors of the 
iz^^ pancake. The average total defect density is defined by = jVAf ^i,z^q{h ^) where 
nq{i, z) = 1 — 5g(i,z),o- (Note that this definition includes defects of both signs.) Since tightly 
bound neutral pairs of disclinations, such as (a) and (c), do not influence local hexatic order 
(i. e. the degree of local six- fold symmetry), it is also useful to define an isolated charge 
density by q*{i,z) = q{i,z) + Y!jq{i,z), where the prime indicates that j runs over the 
in-plane neighbors of pancake {i,z). This definition eliminates bound-pair defects such as 
(a) and (c). The average isolated charge density is then defined by = jq^Hi,z'Kii'^i 
where n*^{i,z) = (1 - 5g*(i,2),o)- 

Once the charges are identified on each plane, the topological defect configuration asso- 
ciated with a given vortex arrangement (lower left in Fig. ^ may be represented as "neutral 
plasma" of line charges of variable lengths Id and both signs (cf. Fig. R|, lower right). To 
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quantify the properties of this plasma, in our simulation, we monitor the distribution of 
defect lengths of either sign V{ld), defined as 

V{l,) = ——{j2Y.l^-n,{t,Zom-n,{z,Zo + l, + l)] U n,{z,z)) (15) 

^^z-lyv q z=zo + l 

We also monitor the in-plane hexatic order parameter ^>q, defined by 

*.~5;;^i:<exp(«%w)> (16) 

where 9ij{z) is angle made by the bond between vortices i and j and the x— axis. The 
results are shown in Fig. ^. Evidently, melting is accompanied by a proliferation of isolated 
dislocation and disclination lines, as well as by a dramatic drop in the in-plane hexatic 
order parameter. Near T^, but still below it, transient line defects are observed to appear 
and disappear, while at there is an abrupt decrease in hexatic order accompanied by 
an almost discontinuous jump in n*^ (cf. Fig. |^). appears to vanish no more than 1 % 
above T^. While we can not rule out a "hexatic line phase" (with hexatic but no long-range 
crystalline order) within this temperature range, our numerical results are consistent with a 
single melting transition where both hexatic and crystalline order simultaneously disappear. 

As the flux line density increases, the interlayer coupling becomes weaker compared 
to the in-plane interaction. Hence, one might expect a dimensional crossover in the defects 
associated with the melting transition. Indeed, we find such a crossover in the solid phase. In 
Fig. ^, we show the length distribution 'P{ld) for 32 layers at several fields {B = 10, 50, 90 kG) 
and temperatures slightly below TmiT/Tm = 0.997,0.982,0.988) selected with the criterion 
of 6R = 0.15 for consistency. In each case, by monitoring C{r,z) and S{c{,z), we verified 
that the system is solid, but very near melting. As evidence that the lattice anisotropy 
depends on field, we note that for the same value of in-plane rms fluctuations {6R = 0.15), 
It is field dependent: at fields of lOkG, 50kG, and 90kG, it is respectively 0.077, 0.100, 
0.114. 

The topological defects have a dramatic cross-over as a function of field. This is obvious 
from the right hand column of the Figure. At the lowest field (lOkG), there are line defects 
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penetrating all the way through the sample, which occur only for temperatures extremely 
close to melting(T/Tm > 0.99). At higher fields, the defect line segments are much shorter 
and occur at somewhat lower values of T/T^. (The latter fiuctuations set in at a lower 
temperatures because they are short and cost less energy to create.) Also shown in gray 
are the defect length distributions for three temperatures above the melting line, namely 
T /Tm = 1.01, 1.04, 1.06 for B= 10, 50, and 90/^(7 respectively. In this case there is no 3D- 
2D crossover: the defect length distributions in these line liquid states are similar for all 
three fields, following a simple exponential form. [However, at the still higher temperatures 
T /Tm ~ 1.05, 1.07, 1.13 respectively (not shown). It > 0.2 and the lines begin to break up 



into pancakes, as cutting and reconnections set in (see Section [1IID| ).] 

To account qualitatively for these results, we note that the defects shown in Fig. |^ 
are line segments of isolated disclinations of various lengths. The energy to create an 
isolated "pancake" (i. e. 2D) disclination defect, as is well known, is proportional to the 
logarithm of the system area. It may be written approximately as Ec ~ Jln(L^/a^), where 
L is the system edge, is the lattice constant of the 2D lattice, and J is some appropriate 
energy, which is of the order of the intervortex interaction energy. The energy to create a 
line of i such disclinations is therefore (neglecting momentarily the interactions between the 
individual pancakes) of the order of E{i) ^ iEc. Since the total number of such defects in 
thermal equilibrium should be proportional to exp{—E{i)/kBT), this argument will give the 
kind of exponential distribution seen numerically in both the solid and liquid phase. 

To further refine this argument, we first note that l/a^ is proportional to the magnetic 
induction B. Hence, for a given area and fixed J, the energy to create a 2D disclination 
increases as InB, which implies that the slope of the exponential dependence should become 
steeper as B increases. This increase is indeed observed in the solid phase, but it seems to be 
much more abrupt than the gradual increase suggested by this argument. Presumably, the 
abruptness stems from interactions between the disclination pancakes. If this interaction 
energy increases exactly as £, the exponential form would be unaltered. Our numerical 
results show some deviation from strictly exponential behavior at low fields, suggesting 
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that the interaction energy is more comphcated. We speculate that this interaction energy 
is stronger at low fields but is reduced at larger fields, possibly by screening from other 
disclinations, leading to the rather sudden transition to short defects seen at high fields. 
But a quantitative theory of this transition remains to be developed. 

Why is there no such transition in the liquid state? A possible but speculative expla- 
nation is that, in the liquid, there are many disclinations (of the order of 0.5 per plaquette 
per layer of the vortex lattice). Hence, they strongly screen one another, and the effective 
sample area in the expression for the creation energy should be replaced by the area 
of a plaquette, which is of order a^. Then the energy to create a line of ^ disclinations is 
independent of field (except possibly through the prefactor J), implying a field-independent 
distribution, as observed numerically. 

Although the above arguments are certainly speculative, the principal numerical re- 
sult - namely, a rather sharp "3D-2D" crossover in the defect structure, may have an ex- 



perimental analog. Obara et al ||32| have recently reported a crossover in multilayers of 
DyBa2Cu307/(Yi_xPrx)Ba2Cu307, in which a 3D vortex lattice showed only 2D correla- 
tions above ~ IQkG. In the experimental sample, of course, the vortex lattice is affected by 
point pins, possibly producing a sharper crossover than we see here. Nonetheless, these pins 
should affect the vortex lattice quite differently at low and high fields, because the defect 
lines are clearly much less rigid at high fields [Q. Specifically, because of this low rigidity 
[as shown in P(/^)], the point pins may cause the defect lines to break into short segments 
rather abruptly at a well-defined magnetic field which may be speculatively identified with 
the transition observed by 



C. Distribution of Magnetic Induction in Solid and Liquid. 

Both the static and dynamic magnetic field distribution are often probed experimentally 
||7|j3^-p6|, using techniques such as ^SR and NMR. To calculate this distribution, we have 



evaluated the instantaneous local magnetic field B{r,t) at a given time t (either by MC or 
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LD), using Clem's prescription ||3^ for the field of a stack of vortex pancakes. In the present 



case, we have an additional complication due to the periodic boundary conditions. This 
difficulty is again solved by including the effects of image pancakes both within the plane 
and along the c-axis. 

To obtain the dynamic evolution of the field distribution, using LD, we consider the 

instantaneous distribution function 

ViB,t) = ^ j dV5{B{v,t)-B), (17) 

which will in general depend on time, as the system approaches equilibrium. For a measure- 
ment which probes the local field averaged over a time t, we obtain the field distribution 
function Vt{B) from the following definition: 

rt{B) = ^JdV6{{B{r,t')),-B), (18) 

where {■ ■ ■)t denotes an average over a time t (either over real time, for a LD simulation, or 
Monte Carlo time). From the limit t — oo, we obtain the static field distribution function 
appropriate for /iSR, (in a typical /iSR, the muons sample the B-field at random points 
in a sample averaged over a typical muon lifetime of ~ 10~^ sec.) using either MC or 
LD. A similar technique has been recently applied to obtain the static field distribution 
in BiSr2Ca2Cu208 using MC |^8[. By varying the duration t over which the local field is 
accumulated, we also find out how rapidly the field distribution approaches the static limit. 

Fig. ^ shows Vt{B) as calculated using LD for t/to ranging from 10^ to 5 x 10^ and three 
different temperatures above and below melting. In the long time limit, the solid phase 
exhibits the characteristically asymmetric profile arising from the static triangular lattice. 
In the liquid, the profile becomes nearly symmetric (and very narrow) because vortices move 
around, producing the same time-averaged field everywhere. The distribution in the solid 
phase is qualitatively similar to that detected by /iSR experiments in BiSr2Ca2Cu208 [0- 
Our symmetric result in the liquid appears to disagree with the experiment. The discrepancy 
may result from the fact that our simulation sample, in contrast to the experimental one, 
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lacks a boundary. The possible influence of such boundaries on the field distribution observed 



by fiSR experiments has been noted and discussed by Schneider et al ||3^. For short time 
scales (t/to < 10^), our distribution retains significant asymmetry even in the liquid phase, 
suggesting that the time-scale for field relaxation is of this order in the liquid. 

Fig. 1^ shows the change in rms width of the field distribution across the melting point. 
In the long-time limit, the temperature dependence of the width as obtained from both MC 



and LD agrees qualitatively with NMR linewidth measurements across melting This 
agreement suggests that the time scale of the NMR measurement is longer than 10^ — lO^to- 
Since the NMR experiments are carried out using a.c. fields of frequency ~ MHz, this implies 
10% < 1/isec or to < lO^i^sec. [This agrees with another estimate of to given below.] Note 
also that our interpretation of NMR linewidth neglects the possibility that the a.c. NMR 
field actually exerts a force on the vortex lattice. 

Two more pieces of information may be extracted from these results. First, the inset in 
the Fig. shows that a MC result taken over 4 x 10^ MC steps gives an rms width closely 
matching the LD result obtained with t/to = 10^. This suggests that a single MC time 
step using our version of the Metropolis algorithm with individual step size Ax = 0^/32 is 
equivalent to ~ 2.5to of LD time. Also, from the LD results with variable t/to^ we can infer 
that melting has a pronounced effect on the field distribution only when the measurement is 
made on time scales longer than ~ lO^to- Iii fact, the required time scale to distinguish solid 
from liquid may be somewhat longer than even this value in the thermodynamic limit. For 
the LD cell size used in these magnetic field calculations, the individual pancakes may have 
5R{t) oc t° with a ~ 1/3 rather than a ~ 1/4 as expected for very long lines (see below). 
Thus, for these thicker samples, we expect that the relevant time scale for field relaxation 
should be closer to (10^)'^/^to ~ 0.5/xsec. 
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D. Relevance of Vortex Line Cutting and Reconnection 



Before turning to further dynamical results, we first consider the validity of neglecting 
vortex line cutting and reconnection. In an earlier numerical study of BiSr2Ca2Cu208 using 



a similar model |11], it was found that including line-breaking effects had little effect on the 
calculated melting properties. In simulations of YBa2Cu307_5 , which has far stiffer lines, it 
is reasonable to expect that the approximation is even sounder. Indeed, even in the liquid 
phase up to T/Tc2{B) ~ 0.97, our numerical results show that the transverse "wandering 
length" It, measured in units of the intervortex separation, is no larger than 0.2. 

To check this approximation another way, define a "vortex collision length" (z by 
^riCz/dY ~ 1- Presumably 5 = 1 in the dilute (low-field) limit, where the transverse 
wandering of a line directed along the c-axis is a random walk with step size It- Using 0.2 
for It, we obtain (z/d ~ 25, which exceeds the thickness of our sample. In the dense regime, 
6 may be smaller than 1 because of the restrictive effects of the repulsive interactions among 
flux lines, leading to an even larger (z/d. 

Cutting and reconnection should occur massively only when the collision length becomes 
comparable to the interlayer spacing, leading to frequent "collisions" of vortex lines. By 
balancing the entropic gain from permutations of vortex connections against the accom- 
panying cost in interlayer coupling energy, we estimate that this condition should be met 
only when It > 0.7, which occurs only well above melting. To verify this, we did two MC 
runs, in one of which we allowed cutting and reconnection according to a Boltzmann weight 
factor obtained from the change in interlayer coupling energy that would be produced upon 
cutting and reconnection. The result shows that this cutting occurs at a negligible rate until 
T/Tm > 1.05, at which point about 12% of recombination attempts are accepted. It at this 
temperature was found to be about 0.16, a value which may roughly be taken as a kind of 
"Lindemann melting criterion" for flux cutting in the liquid state. Thus, over much of the 
liquid regime, we conclude that the thermodynamics of this model can be treated without 
considering flux line breaking and reconnection. This conclusion justifies our treatment of 
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flux line dynamics in the same approximation. 



E. Langevin Dynamics of Vortex Line Liquid and Solid; Slow and Fast Relaxation 

Having justified our non-breaking assumption, we return to the dynamics of vortex hne 
liquids and solids. We begin by considering the motion of a single vortex line within the 
overdamped dynamics of eq. (P). In the limit of a very long line in which each pancake is 
subject only to thermal Langevin noise and a harmonic interlayer restoring force, eq. (P) can 
be solved analytically. The derivation is similar to that of Ref. P] for more general cases. 
The result for the mean square displacement of a pancake from its initial position at time t 
is 

which predicts a su6-/mear time-dependence By contrast, an uncoupled pancake vortex 
has an ordinary diffusive transverse motion, in which (|rj^fc(t) — ^(O)^) oc t. 

To test this behavior, we have examined the long-time behavior of the quantity 6R{t) 
defined in eq. ([TT|) from LD in the limits of (i) a single pancake vortex; (ii) a single long line; 
(iii) single short line segments; and (iv) an ensemble of 64 lines in 8 and 16 layers, at various 
temperatures both below and above T^. For the single pancake, we observe ordinary diffusive 
behavior. From the limit t/to —>■ 10^^, we deduce T> = lAx 10~'^cm2/sec for T = QOK (using 
to = 1.2 ■ 10~^^sec; see below). This diffusion constant seems to agree reasonably well with 



experimental estimates for an extremely anisotropic system such as BiSr2Ca2Cu208 |^0|- ^ 
single line of up to 1000 vortices shows SR{t) oc t^^^ as predicted analytically, while lines 
shorter than Ad show behaviors close to a 2D diffusion. 

Fig. p shows case (iv) for B = 50kG and 16 layers of YBa2Cu307_5 at five temperatures 
above and below T^. For T < T^, vortex excursions are limited to radii smaller than 
the appropriate Lindemann distance of around 0.20^. For T > Tm, the rms displacement 
seems to grow with an approximately t^^^ behavior as expected from the analytic estimate 
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for a single line. In the case of an eight-layer system(not shown), we observed SR{t) oc 
t^/^, possibly indicating that the system lies between the long-line and 2D limits. In both 
cases, the motion of individual pancakes is slower than in the usual Brownian diffusion. 
Hence, in the long time limit, a diffusion constant defined on the assumption of a linear 
t-dependence will vanish even in the liquid regime and even without pinning. At sufficiently 
high temperatures, of course, line cutting will eventually set in (though not in the present 
approximate model). With line cutting, the system may then cross over to a 2D liquid with 
ordinary diffusive behavior. Therefore, there is an interesting possibility of two different 
types of liquids characterized by different diffusive behavior, with a smooth crossover between 
them. 

As a further means of studying relaxation in the solid and liquid phases, we have mon- 
itored the LD evolution of various thermodynamic quantities, such as various components 
of internal energy and the defect density, after the system is initialized in some arbitrary 
non-equilibrium state. For each temperature, we considered five different initial states. Each 
initial state was prepared by randomly and rigidly displacing the individual vortex lines from 
a perfect triangular lattice by an amount not exceeding a fraction of the mean intervortex 
spacing. In most rapid exponential relaxation of energy is observed. By fitting the 

time-dependent internal energy to the form a + 6exp[— t/r^.], we find a relaxation time 
which increases with increasing temperature (cf. Fig. ^ for which = 0.31). This increase 
is once again due mainly to an energy scale which diminishes [or a Xab{T, B) which increases] 
with increasing T. 

Superimposed on this overall trend, there may be a peak in corresponding to enhanced 
fluctuations precisely at the melting temperature (although there are large uncertainties at 
this temperature). This behavior is consistent with expectations for a phase transition with 
a kinetic or "critical" slowing down, but the feature is generally concealed by the steep 
increase of with increasing temperatures, and can hardly be used to distinguish between 
possible first-order and continuous phase transitions. The middle row of Fig. |^ shows the 
progression of the density correlation function from an initial random configuration to a well 
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ordered state over a period of about lOOto at a temperature T /Tm = 0.9. The lower row 
shows the local instantaneous magnetic field profile for the initial and final configurations 
in a particular layer. The cores of the vortices lie at the centers of the bright regions. 

In Fig. [ly, we show a typical relaxation of total internal energy over time after an initial 
randomization, this time with = 0.63. Also shown is the density of isolated disclinations, 
n*^. In contrast to the internal energy, which relaxes exponentially, the disclination density 
relaxes toward zero (i. e., an elastic lattice) with a roughly hit behavior. The disclinations 
finally disappear (via pair annihilation) long after the internal energy has nearly equilibrated. 
The spikes in n*^ on top of the overall logarithmic decay may correspond to a large-scale 
rearrangements of vortices, possibly involving activated processes. The energy relaxation 
itself shows no obvious signature of any such activated processes. 

To understand real materials, we next examine how point pins infiuence these general 



features. In Fig. |TT], we show the relaxation of the in-plane component of the internal energy 
for different pinning strengths (0 < < 5) with an areal density equivalent to Bp = 888 
kG. The relaxation is qualitatively different from the weak-pinning cases (a^ < 3), in which 
the energy typically decays exponentially as in clean systems. At ap = 5, for example, the 
in-plane part of the energy varies logarithmically with time almost from t = 0. 

For pins of any strength, we find that the system never relaxes back to the perfect 
triangular lattice within our simulation times; a significant fraction of disclinations always 



survives even after t/to = 5 x 10^ (cf. lower part of Fig. [T^). Possibly, early-stage relaxation 



is controlled mainly by pinning forces. But once the pinning energy is nearly optimized by 
this quick relaxation, all driving forces for relaxation become comparable in strength (after 
t/to ^ 3 X 10^), and relaxation slows dramatically. Any further relaxation must reduce 
the elastic strain energy without sacrificing pinning energy. To accomplish this, both total 
and isolated defect densities slowly but steadily decrease (with many fiuctuations). Also, 
n^/rid [lower panel of Figure] decreases noticeably for t/to > 3 x 10^. In other words, 
energy relaxation in this regime is proceeding via a decrease in the number of isolated defect 
configurations. 
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The slow elastic relaxation just described must involve vortex rearrangement. It may 
therefore affect measurements of the local magnetic field and the total magnetization. In- 
deed, such quantities have long been known to depend strongly on the time scale over which 
they are measured. Thus, our model can help to explain such slow relaxation of magneti- 
zation, in disordered samples. The calculations also suggest that the slow relaxation results 
predominantly from the rearrangement and annealing-out of topological defects. 

Next, we briefly discuss the value of to, which is crucial in connecting our numer- 
ical results to experiments on real materials. The friction coefficient rj is related to 
the flux flow resistivity p// by the equation pjf = Assuming a flux- flow resistivity 
P//(B = 50kG, T/Tm = 0.9) ~ 0.16 p„, where p„ is the normal state resistivity, and estimat- 
ing pn as its value at Tc(0), p„(Tc(0)) ~ lOOpfi-cm |^ for a single crystal of YBa2Cu307_5 , 
we obtain pff ^ 1.6 x 10"^'' sec. The corresponding value of to, assuming the time-step in 
Sec. lie, is to ~ 10^^'^ — 10~^^ sec. Since our value of p// is probably an overestimate, this 



value should be taken as the lower bound for to ||42[] . 

Finally, we briefly return to the influence of pins on the dynamics of flux lines, in light of 
this value for to. In real materials, such pins will strongly affect flux line dynamics, typically 
slowing down their motion by several orders of magnitudes (as we have noted above). This 
has some implications for interpreting the experimental results which are sensitive to the 
time evolution of local magnetic field. One such experiment is the spin-echo NMR probe 
40| , |43[| in the highly anisotropic Tl-based compounds. This has been interpreted as implying 



a vortex diffusion constant T> of 10^^ ~ 10~^cm^/sec at temperatures T ~ T^, and of ~ 10^^ 
- 10~^cm^/sec for T > T^. On the other hand, more recent spin-echo results on aligned 
powders of YBa2Cu307_5 seem to show lack of diffusion even on very long time scale |]35| , |36| . 
These results may indicate the presence of strong pins in the powder samples, which greatly 
slow down the dynamics. Another possibility is that the discrepancy is somehow due to the 
nondiffusive behavior of long lines discussed above. 
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IV. DISCUSSION AND CONCLUSIONS. 



We have presented a detailed numerical study of the vortex structure in a model for 
YBa2Cu307^5 , both above and below the flux lattice melting temperature, using Monte 
Carlo and Langevin simulations. Our results indicate that in clean YBa2Cu307_5 , there is 
a weakly first-order melting transition, with a very small heat of fusion which agrees the 
calculations of at similar fields. 

We also find that there is a striking change in distribution of local magnetic fields at 
melting: the time-averaged magnetic field sensed by a fixed particle in the sample has a 
very narrow distribution in the liquid (because the vortices move around in the liquid), but 
acquires an asymmetric distribution of finite and temperature-dependent width in the solid 
phase. This behavior agrees with both /iSR and NMR experiments in YBa2Cu307_,5 . Al- 
though some previous MC calculations have suggested a similar transition, our LD sim- 
ulations provide more information. In particular, because the time constant of the Langevin 
simulations can be extracted from measurements of the liquid-state resistivity, we can esti- 
mate the measurement time necessary to distinguish the solid and liquid field distributions. 
This time appears to be about 0.5 /isec in YBa2Cu307^5 at a field of 50 kG. 

Our calculations show that there is a qualitative change in the structure of the topologi- 
cal defects in the crystalline phase just below melting. Namely, even in clean YBa2Cu307_5 , 
there are long lines of disclinations parallel to the c axis at relatively low fields (B ~ lOkG), 
but these break up into short ("2D") disclinations at higher fields. We believe that this 
breakup may be an essential ingredient in the "3D-2D" transition recently reported in 
BSCCO [^,^,0]. In clean systems, we find numerically that the distribution of defect 
line lengths is approximately exponential in the liquid, but deviates from exponential in the 
solid phase. We account for these forms by a simple model of topological defects in the 
vortex lattice of a layered superconductor. 

A remarkable result of our Langevin simulations is that the vortex lines move nondiffu- 
sively even in the liquid state - that is, the rms displacements of the vortex pancakes vary 
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like t", where a < 0.5. This agrees with the exact result derived for a single long line |P9|| , 
that the rms displacement of a sufficiently long vortex line subject to thermal noise varies 
as t^^^. Our simulations suggest that the behavior is also found in a dense line liquid with 
strong short range repulsive interactions. While the full transport implications remain to 
be explored, this result is consistent with recent NMR measurements on YBa2Cu307_5 PB 
which also suggest that in the liquid state vortex lines move sub diffusively. 

The LD calculations also reveal that the vortices relax much more slowly as T ^ Tc2{B). 
Such slow relaxation is actually built into the TDGL equation, whose time constant is 
proportional to the deviation of the Ginzburg-Landau free energy from its minimum value. 
Since the minimum becomes ever shallower as Tc2{B) is approached, the relaxation time 
constant should become ever longer. It would be most interesting to have experimental 
confirmation of this behavior. The Langevin results also suggest that the time constant may 
increase as the (first-order) melting transition is approached from either side. 

A final conclusion has to do with the infiuence of point pins on the dynamics. Our 
Langevin results show that, after any disturbance, the energy relaxes nearly exponentially 
back to equilibrium in a pin-free sample. With strong pins present, however, the energy 
relaxation is close to logarithmic in time. At even longer time scales, there is a crossover to 
a logarithmic relaxation with a different slope. We identify this slower relaxation with the 
logarithmically slow annealing out of topological defects in the disordered samples (these 
same defects are also slow to disappear in the pin-free case). It has long been known 
that many properties of the high-T^ materials vary logarithmically slowly with time (most 
notably, the sample magnetization). Our model, which treats in a fairly realistic manner 
a disordered YBa2Cu307_5 sample, clearly produces such slow relaxation, and identifies it 
with a particular type of topological defect dynamics. 

In summary, we have presented a detailed study of the nature of the fiux line melting 
based on the vortex representation of the Lawrence- Doniach model, using both Monte Carlo 
and Langevin simulations in a consistent way. We calculate a wide range of properties of 
both the vortex solid and liquid phases which are consistent with experiment, and which 
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shed new light on the topological defects which underhe the melting process, as well as the 
time-dependent magnetization and magnetic field distribution, of these materials. 
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FIGURES 

FIG. 1. Probability distribution V{U) for the total internal energy Utotai{T) per pancake, av- 
eraged over two different time windows (a) and (b) (specified in Fig. |2|), at the melting temperature 
Tm- Calculation is carried out for YBa2Cu307_5 at B = 50kG. Uhex is the corresponding energy 
for a perfect triangular lattice of straight vortex lines (the minimum energy configuration). The 
distributions in each window are obtained by dividing the data from the Monte Carlo configura- 
tions shown in Fig. Q into 256 energy bins, each of width MJ/ksT = 1.4 x 10"^. The curves are 
least-squares fits of these data to Gaussian distributions. The error bars at various values of Utotai 
are the rms deviations of the values of V{U) inferred from ten bins in the vicinity of U . Inset: 
in-plane rms displacement c^, of pancake vortices from their equilibrium lattice positions in units 
of lattice spacing, as calculated using both Monte Carlo and Langevin techniques. 

FIG. 2. Evolution of the in-plane component of internal energy, Uinpiane{T), at the melting 
temperature, for B = 50kG as calculated via MC simulation. Insets: Density-density correlation 
C(r,0), averaged over the two different time windows (b) and (c). These plots show that the 
system is alternating between lattice and liquid phases. 

FIG. 3. Snapshot of disclination distribution in a single layer of a vortex liquid containing 
256 lines, as determined by Delaunay triangulation. Black dots: fivefold disclinations; gray dots: 
sevenfold disclinations. (a) A pair of bound disclinations, equivalent to a dislocation; (b) an isolated 
disclination; (c) a pair of bound dislocations. The bottom row shows a typical vortex line liquid 
and its representation in terms of disclinations of either sign. 

FIG. 4. Hexatic order parameter |V'6p (left scale) and density of isolated disclinations (right 
scale) as a function of temperature in a 64-line system of thickness 8 layers with B = 50kG. 
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FIG. 5. Darker symbols: probability distribution of disclinations of length at B = 10, 50, 
and 90 kG for a 32-layer system, at temperatures such that 6R{B,T) = 0.15, i. e., slightly below 
melting. Lighter symbols: corresponding probability distributions for the same fields and system 
at temperatures deep into the liquid phase. In the right column, snapshots of typical defect 
configurations in the solid phase are shown for each field. Black and gray represent disclinations 
of either sign. 

FIG. 6. Dependence of time-averaged magnetic field distribution Vt{B) on time window t/tQ 
used for average, at an applied field of 50A;G. The time evolution is calculated at several temper- 
atures from Langevin dynamics using Clem's prescription for computing the magnetic field ||37|| . 
The distribution is plotted as a function of B/Bav — 1, where Bav is the space- aver aged magnetic 
induction. Successive distributions in each vertical panel are displaced horizontally by 0.02 units. 

FIG. 7. Mean-square width {dB^/Bl^ of the magnetic field distribution V{B/Bav), plotted 
as a function of T/Tm{B) at an applied field B = 50A;G. The width is calculated from Langevin 
dynamics for time intervals of varying durations: t/to = 10^,10^,10^ (inset), and 5 x 10^ (main 
Figure). Also shown in the inset is the mean-square width as obtained from a MC simulation 
sampled over 4 x 10-^ MC steps. 

FIG. 8. (a) Rms transverse displacement 6R{t) of vortex pancakes at several temperatures 
both above and below the melting temperature T^, as calculated by Langevin dynamics at B = 
50kG for Nz = 16, and plotted in units of a^, the mean vortex separation. For T < T^, the rms 
displacement saturates at a value approaching the Lindemann number cl ~ 0.18 at T^. In the 
liquid regime, the displacement increases roughly as t^^^, as expected for a long line, and not as 
t^^'^ as expected of diffusing particles, (b) Root-mean-square "wandering length" It at the same 
temperatures. Note that It saturates even in the liquid state, though its value continues to increase, 
primarily because of the reduction of line tension with increasing temperature. 
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FIG. 9. Top: relaxation time Tr describing return of vortex system to its equilibrium configu- 
ration, following an initial perturbation, as plotted for several temperatures in the flux solid and 
flux liquid state. Calculation is carried out as described in the text. Error bars are standard de- 
viations from five different initial configurations. Middle: snapshots of density-density correlation 
function C{r,z = 0) for three representative configurations during relaxation at T/Tjn = 0.89. 
Bottom: evolution of local field B(r,t) in a specific ab plane of the sample, for the first and third 
configurations of middle panel. 

FIG. 10. Time-dependence of total internal energy Utotal{T^) per pancake (left-hand scale) 
and the isolated disclination density (right-hand scale), as calculated via LD simulation at 
B = 50kG, T/Tm = 0.898, starting from a randomized initial configuration in a pin-free sample. 
U/ieo; is the internal energy of an ideal hexagonal lattice at the same temperature. 

FIG. 11. Relaxation of in-plane portion of internal energy per pancake, \Jinpiane, for point 
pins of varying strength(0 < ftp < 5). The initial configuration is a randomized distribution of 
straight vortex lines with = 0.63 (see text). Note the gradual crossover from an exponential 
to a logarithmic time dependence as the pinning strength increases. At the largest value of ap, 
the relaxation is dominated by changes in the pinning energy. Other conditions same as in the 
previous Figure. Uhex is here the in-plane part of the internal energy of an ideal hexagonal lattice 
(per pancake) at the same temperature. 

FIG. 12. Time-dependence of interlayer part of internal energy \Jj(T), pinning energy Upi„(T), 
in-plane component of elastic energy ^inpianei'^)^ ratio of isolated disclination density to 
total disclination density n^, for strong pinning {ap = 5), at a temperature of T/Tj^ = 0.898 where 
Tm is the melting temperature for a pin-free sample. 
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